function h = plot_q(grid,base,Q,varargin)
% h = plot_q(grid,base,Q)
%
% Plot a vector field Q.
%
% h = plot_q(grid,base,Q,scale_factor)
% Scale the arrows with scale_factor
%
% h = plot_q(grid,base,Q,scale_factor,0)
% Plot arrows only at the barycenters

 if(nargin>3)
   scale_factor = varargin{1};
 else
   scale_factor = 1; % no scaling
 end

 if(nargin>4)
   meshg.pg = 1/(grid.d+1)*ones(grid.d,1);
 else
   meshg = load('plot_mesh.mat');
 end

 % evaluate the basis on the meshg.pg nodes
 OMEGA = zeros(base.me.d,base.mk,size(meshg.pg,2));
 for i = 1:base.mk
   for j = 1:base.me.d
     OMEGA(j,i,:) = permute(ev_pol(base.o_s{j,i},meshg.pg),[3 1 2]);
   end
 end

 Xp = zeros(grid.d,grid.ne,size(meshg.pg,2));
 Qp = zeros(grid.d,grid.ne,size(meshg.pg,2));
 
 for ie = 1:grid.ne
    
   X_plot = grid.e(ie).b * meshg.pg;
   for j = 1:grid.d
     X_plot(j,:) = X_plot(j,:) + grid.e(ie).x0(j);
   end

   Q_el = Q((ie-1)*base.mk+1 : ie*base.mk);
   for j = 1:grid.d
     Q_plot(j,:) = Q_el' * permute(OMEGA(j,:,:),[2 3 1]);
   end

   for j = 1:grid.d
     Xp(j,ie,:) = X_plot(j,:);
   end
   Q_plot = 1/(2*grid.e(ie).vol)*grid.e(ie).b*Q_plot;
   for j = 1:grid.d
     Qp(j,ie,:) = scale_factor*Q_plot(j,:);
   end
   
 end

 h = figure;
 switch grid.d
  case 2
   pdemesh(h,grid,'k');
  case 3
   pdemesh(h,grid,'y');
  otherwise
   error(['Dimension ',num2str(grid.d),' not supported']);
 endswitch
 hold on
 
 switch grid.d
  case 2
   for ie = 1:grid.ne
     hh = quiver(Xp(1,ie,:),Xp(2,ie,:),Qp(1,ie,:),Qp(2,ie,:),0);
     set(hh,'color',[rand(1),rand(1),rand(1)]);
   end
  case 3
   for ie = 1:grid.ne
     hh = quiver3(Xp(1,ie,:),Xp(2,ie,:),Xp(3,ie,:), ...
                  Qp(1,ie,:),Qp(2,ie,:),Qp(3,ie,:), 0);
     set(hh,'color',[rand(1),rand(1),rand(1)]);
   end
 endswitch

 title('q_h');
 axis equal

 hold off

return

